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The residual deposits usually left near the contact line after pinned sessile colloidal droplet evaporation are 
commonly known as a "coft^ee-ring" effect. However, there were scarce attempts to simulate the effect, and 
the realistic fully three-dimensional (3D) model is lacking since the complex drying process seems to limit 
the further investigation. Here we develop a stochastic method to model the particle deposition in 
evaporating a pinned sessile colloidal droplet. The 3D Monte Carlo model is developed in the 
spherical-cap-shaped droplet. In the algorithm, the analytical equations of fluid flow are used to calculate the 
probability distributions for the biased random walk, associated with the drift-diffusion equations. We 
obtain the 3D coffee-ring structures as the final results of the simulation and analyze the dependence of the 
ring profile on the particle volumetric concentration and sticking probability. 

The coffee-ring effect has been observed in various experiments with different colloidal suspensions^"^^, 
including polymers^"^, quantum dots^, nanoparticles^"^^, bacteria^ \ and biofluids^^"^^, such as urine^^, 
blood^^ and serum^^. The effect has initially been investigated by Deegan et al/^ and explained as a 
ubiquitous phenomenon caused by the replenishing outward capillary flow and uneven evaporation rate^^'^^. 
The previous attempts to simulate the effect were commonly based on the numerical techniques, such as the 
finite-element method^°, applying the analytical equations to figure out the evolution of the average particle 
density profile. Another family of approaches, based on the Monte Carlo methods, uses randomly distributed 
particles on the discrete lattice and calculates the probabilities of the possible particle motion directions on each 
simulation step. One recent attempt to build a Monte Carlo model of the droplet evaporation, capillary flow and 
contact line deposition has been performed by Kim et al.^\ The model used a simple power-law assumption about 
the particle motion on the basis of the flow analysis of Deegan et al.^^ and reproduced the contact line deposition 
profile^. However, the simulation has been limited to the simplified radial particle motion and has not considered 
the vertical velocity component and dependence on the vertical coordinate. In the work of Yunker et al.^^, the 
coffee-ring growth was modeled in two dimensions (2D) as the Poisson-like process of random particle depos- 
ition coupled with the surface diffusion. The model, focusing only on the vicinity of the contact line, has not 
considered the droplet evaporation and inward flow dynamics. The Monte Carlo approach based on the biased 
random walk (BRW) has been recently used to investigate the transition from the coffee-ring deposition towards 
the uniform coverage in 2D^'^'^^. In the present work we further develop the BRW method into the more realistic 
3D domain and use the flow analysis of Hu and Larsson^^"^^ to calculate the corresponding probabilities of the 
sampled particle moves on each Monte Carlo Step (MCS). This advancement allows achieving a full 3D structure 
of the coffee ring and analyzing the thickness profile of the structure. The evolution of the ring shape and 
dimensions is observed during the entire time period of the droplet drying. 

The mathematical model is based on the BRW approach^^"^^ in a 3D domain with a cubical lattice structure, 
which mimics the shape of a pinned sessile droplet (see Fig. 1) with a spherical cap with a pre-defmed initial 
contact angle value, Oq, and apex height, hg. The domain boundary surface continuously shrinks in the vertical 
direction during evaporation. The geometrical shape is assumed to remain spherical cap, and the contact angle 
value continuously decreases with time. It is assumed for simplicity that the droplet apex height reduces linearly 
with dimensionless time, P^, /z(Oj) = ho{l — t). At the initial state, the particles are uniformly distributed inside 
the domain according to the defined volumetric concentration, (j). Afterwards, on each MCS of the simulation, 
each particle may perform a biased random move into a vacant cell within the domain. Each particle move is 
sampled from the seven- component finite set, where six elements correspond to the six possible directions of 
motion, and the seventh possible move is "waiting", or keeping particle position in the lattice until the next MCS. 
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Figure 1 | The schematic representation of the 3D model 
implementation. The spherical cap-like volume of the cubical lattice is 
filled with liquid (blue) and randomly distributed particles (red cubes), 
(a) shows the isometric view of the modeled droplet, and (b) shows the side 
view of the droplet. 

To calculate the probabilities corresponding to the chance of drawing 
a given direction for a particle, first of all, the dimensionless velocity 
components Vx, Vy, Vz are calculated from the analytical expres- 
sions^^'^^ for each particle (see the details in the Methods section). 
The velocity values significantly depend on time and the particle 
position, rising near the three-phase domain boundary. Mean- 
while, the lattice configuration restricts the simulated particle dis- 
placement with one lattice cell at one MCS. Therefore, the duration of 
each MCS might be different, and the tick size of the clock for t is 
adjusted to the fastest moving particle. The drift velocities of the 
slower particles are scaled by the factor of velocity of the fastest 
particle. Thus all the rescaled velocities fall in the range between 0 
and 1, and the preliminary probability distribution of the sampled 
move directions can be calculated. However, one more important 
modification implementing the "waiting" moves is required to keep 
the diffusion parameter consistent in a BRW model with a strong 
bias^^'^^. It has to be remembered that the random Brownian (or un- 
biased) displacement is dependent on the duration of MCS. 
Therefore, shorter simulation steps corresponding to the faster drift 
velocities should have lower probabilities of moving in any direction, 
different from the direction of v. Lowering the diffusion component 
without changing the drift component results in increasing the wait- 
ing probability during the current MCS^^. Finally, seven probability 
components are re-calculated for each particle on each MCS, 
Pmove = {px+ ,px- ,py+ ,py^ ,pz+ ,pz^ ,pw}'^ the subscripts X, )/, z indicate 
the three Cartesian axes, " + " and " signs indicate whether this 
direction matches the direction of v projection on this axis or is 
opposite, andp^ is the probability of the "waiting" move. Oncep^ove 
is determined, the direction of the current move is sampled from the 
biased distribution. 

An additional method is used to improve the convergence of the 
simulation run. Since the velocity value near the three-phase line 
region becomes very high, especially at later stages, the direct scaling 
by the maximum velocity may lead to the MCS duration approaching 
zero, and hence reaching the end of the evaporation process dX ~t=\ 
may require an infinite number of MCS. To avoid this divergence, the 
particles reaching the three-phase contact line or inner side of the 
coffee ring are excluded from calculations of the maximum velocity. 
It is reasonable, since the outwardly moving particles reaching the 



domain boundary are not able to move further, and their actual 
velocity becomes zero. The same situation happens with the particles 
reaching the boundary coffee-ring composed from the blocked part- 
icles at the later stages. Even if the analytical equations show a high v 
value for such a particle, the model rules do not allow it to move into 
the occupied cells, and the particle may be excluded from the simu- 
lated dynamics. Therefore, the duration of current MCS is calculated 
based on the velocities of all the other particles which have not yet 
been blocked by the boundaries. 

The model is also considering the possible particle agglomeration 
into clusters on each MCS. An additional parameter, psuch controls 
the sticking probability of the particles. The mechanics of the process 
is implemented according to the Diffusion Limited Aggregation 
(DLA) rules^'^'^^ The aggregated particles continue their further 
motion as one inseparable object (cluster), and their averaged posi- 
tion is used to calculate pmove and sample the corresponding dir- 
ection. The aggregation event is the binary random variable 
sampled for the particles occupying the neighbouring cells with the 
probability equal to p stick' {^i)^ where p stick is fixed during the run, 
and At is the duration of current MCS. Particles reasonably have 
higher aggregation chance during the longer time steps. 

Results 

The simulations are performed using the following values of the fixed 
parameters: the droplet base diameter is set to 500 cells; the initial 
contact angle value, Oq, is assumed to be 7r/4; the non-dimensional 
particle size is 1 X 1 X 1 cells; and the dimensionless initial evap- 
orative flux at the droplet apex is fixed at 7(0,0) = 1. The volumetric 
particle concentration, (j), and sticking parameter, p^^/^icj are varied in 
consecutive simulation runs to investigate the influence on the cof- 
fee-ring shape and final particle distribution. Figure 2 presents the 
results of a sample simulation run with at 0.01. The image sequence 
shows the progression of the droplet evaporation and the coffee-ring 
forming process. The process starts with the uniform particle distri- 
bution inside the spherical cap (see Fig. 2 (a)), then the droplet 
consequently shrinks while the ring starts to build up (Figs. 2 (b- 
h)); at the later stages the rest of the particles joins the ring structure 
and completes the inner ring side (Fig. 2 (i-r)), finally only very few 
particles are left inside the ring, and the droplet completely dries out 
(Fig. 2 (s-t)). The final image of the ring well resembles some pre- 
vious experimental observations^. 

The investigation of the coffee- ring shape dependence on the par- 
ticle concentration and sticking probability is performed as well. 
Figure 3 displays the final results of the multiple simulations with 
the different particle concentrations varying from 0.01 to 0.05, and 
P stick values at 0.1 (a-b rows in Fig. 3) and 1.0 (c-d rows). Image 
collection (Fig. 3 (al-a5)) presents a gradual increase of the ring 
width with the rising number of particles in case of low sticking p^^^^^ 
= 0.1. Figure 3 (bl-b5) shows that only small clusters, indicated by 
different color, are formed in these cases, and most of the particles do 
not stick together. The shape of the ring looks generally smooth, and 
the width rapidly grows with increasing (j). However, increasing the 
sticking chances ten-fold to p stick ^ 1 provides visible changes in the 
final ring structure for higher particle concentrations (Fig. 3 (c3-c5)). 
More particles stick together to form slower moving clusters, thus 
slightly distorting the inner ring boundary and increasing the particle 
presence near the domain center. Cluster diagrams in Fig. 3 (dl-d5) 
show that p stick = 1 is high enough to cause forming of a large ring 
aggregate for (j) > 0.01. Differently, the coffee rings in insets Fig. 3 
(al-a5) for p stick = 0.1 are composed from multiple small clusters or 
separate particles (Fig. 3 (bl-b5)). The diagrams in Fig. 4 present 
thickness profiles of the final results corresponding to Fig. 3. The 
comparison between color maps shows that thickness of ring struc- 
tures is distributed rather similarly for p stick = 0.1 (shown in Fig. 4 
(al-a5)) and pstick = 1 (shown in Fig. 4 (bl-b5)). However, the 
thickness of the rings with a low sticking changes very gradually. 
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Figure 2 | The progression of the single droplet drying process over time in a Monte Carlo simulation run. (a) shows the initial system state a.tt = 0.001 
from isometric view (top image) and side view (bottom image); (b-t) show the corresponding images of the drying droplet after (b) ^ = 0.1, 
(c) t = 0.2, (d) t = 0.3, (e) t = OA, (f) t = 0.5, (g) t = 0.6, (h) 1 = 0.7, (i) 1 = 0.8, (j) 1 = 0.85, (k) 1 = 0.9, (1) ^ = 0.91, (m) ^ = 0.92, (n) 1 = 0.93, (o) 1 = 0.94, 
(p) ^ = 0.96, (q) 1 = 0.97, (r) 1 = 0.98, (s) ^ = 0.99, (t) 1= 1. The frill animations are available as the Supplementary Video- 1 for the isometric view and 
Supplementary Video -2 for the side view. 



reaching the peak closer to the three-phase line and decreasing 
towards the domain center (Fig. 4 (a4-a5)). Meanwhile, higher stick- 
ing leads to the emergence of more complex branched aggregates for 
(j) > 0.03 (Fig. 4 (b4-b5)), thus slightly distorting the profile of the 
ring. The corresponding curves in Fig. 4 (a6-b6) indicate a higher 
particle presence inside the ring for high-sticking and high-concen- 
tration results (Fig. 4 (b4-b5)). 

Discussion 

The presented algorithm might serve as a simplistic basis for devel- 
oping more specific or sophisticated models, potentially including 
the effects of Marangoni convection, solution viscosity, etc. There are 
several additional important rules considering the particle 
mechanics in the presented model, worth mentioning: 1. The part- 
icles are not allowed to overlap with each other (conservation of the 
global particle concentration, (/>). This rule allows particles to stack 
onto each other and form a solid structure. 2. Once the particle leaves 
the droplet domain as a result of the droplet shrinking, the only 
allowed move is in the vertical down direction, unless the underlying 
cell is occupied by the other particle or substrate. This restriction 
prevents the particle diffusion in the absence of liquid but allows the 



particles drop due to the gravity. 3. All the moves outside the current 
domain are prohibited (the barrier boundary condition). Thus, the 
simulation starts at f = 0 with the particles randomly distributed 
inside the volume of the spherical cap with h=l and finishes at 
t=l, h = 0 with the particles immobilized on the empty substrate 
or atop the other particles. 

To summarize, the coffee-ring effect is simulated in the 3D 
domain, while including both particle diffusion and the capillary 
particle transport towards the three-phase contact line. The spatial 
domain of the model is the evaporating sessile droplet represented as 
a spherical cap. The time scale of the simulation is the full drying 
process, from the placement of the drop onto a substrate till the final 
dry- out of the remaining solvent. The simulation result provides a 
full 3D model of the coffee ring, and the vertical ring thickness profile 
is measured. The final shape of the simulated ring structure shows a 
reasonable dependence on the volumetric particle concentration and 
colloidal aggregation parameter. 

Methods 

The distribution pmove = {px+ ,px^ ,py+ ,py_ ,pz+ ,pz^ ,pw} for each particle is deduced 
from the dimensionless velocity components as 
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Figure 3 | Aggregated simulation results showing the final coffee-ring structures for different volumetric particle concentrations and sticking 
parameter values. (al-a5) show the final dried states corresponding to particle concentrations, 0, at (al) 0.01, (a2) 0.02, (a3) 0.03, (a4) 0.04, (a5) 0.05, and 
sticking parameter value pstick = 0. 1. The remaining particles are shown in red color. (bl-b5) show the cluster diagrams corresponding to the results from 
(al-a5). Different particle clusters are indicated by different colors. The small size of the clusters is caused by a very low particle sticking probability. 
(cl-c5) show the final dried states corresponding to particle concentrations, 0, at (cl) 0.01, (c2) 0.02, (c3) 0.03, (c4) 0.04, (c5) 0.05, and sticking 
parameter value p5ftcit — 1, particles are shown in red. (dl-d5) show the cluster diagrams corresponding to the results from (cl-c5). The main, ring-shaped 
aggregate is seen in (d2-d5) and indicated with green color. Other, smaller clusters are seen inside the main ring and are more evident for the highest 
particle concentration in (d5). 



Py--- 
Py.-- 

Pz--- 
Pw=l-P: 



' 2N 



' 2N 



-Py- 



-Py+ -Pz- -Pz+ 



(1) 



where AT = 3 for 3D space, Vmax > 1 is the dimensionless velocity of the fastest particle, 
Vx, Vy and Vz are the dimensionless velocity components, calculated from the ana- 
lytical equations for the liquid flow^^'^^: 

v.(r,zJ)=iA[(l-F)-(l-F)-^^)](|-2l)-#(|^ 
v,(r,zj) = 1^ + -r^)-^-('^)-^] - f) + f ^ [(1 -F) 

- - 1^)^(0'^) + - f) + - f ) - ^^P(o,^), 

= -2/ior(/(0)i(0)(l + 1) -/(Ma,Ar,v), 



(2) 



where the dimensionless variables are defined: radial velocity Vy = '^r'tf / Ry vertical velocity 

= ^ztf/ Ry time t = t/tf, radial coordinate f = r/R, vertical coordinate z = z/ho, and 
liquid layer thickness h = h/hQ,m which tf is the total drying time, and R is the droplet 
base radius, /i{6) is linearly dependent on the current contact angle value 6, /i{6) =1/2 
— 9/n; g' is the first-order partial derivative of the function, g{f,t), with r; / = J{0,t) is 
the dimensionless evaporative flux at the top of the droplet surface, 

Joe (0.270^ + 1.3) 1^0.6381-0.2239(0- according to Hu and Larsson^^ and 

J{Ma, AT, v) is a Marangoni effect function of the cumulative influence of the 



temperature and surfactant gradients. Our simulations show that low dimensionless 
values off if < 5) do not significantly affect the final structures, therefore /term is 
neglected in the presented results (/= 0). The influence of the strong Marangoni 
effect (/> 5) on the coffee-ring structures in 3D is a subject for future study. The 
expressions of and are calculated for each particle and recalculated for each MCS, 
then the horizontal velocity components in x andy directions are simply derived from 
projections = v^xlr and Vy = Vrylr. 

The maximal absolute value of the particle velocity Vmax on a given MCS is taken as 
maximum max(yx,Vy,Vz) from all the calculated particle velocity components. 
Particles reaching the domain boundary or the coffee ring are excluded from this 
calculation, since their actual velocity vanishes. 

Cluster-cluster aggregation is performed by examining each individual cluster for 
neighbouring particles from the other clusters. If the particles are present in the 
neighbouring cells, the binary variable is sampled with probability oipstick'{Ai), At is 
calculated from At= 1/ {NBVynax)y in which B is the number of cells in the lattice 
radius (250 in the current work). If the sampled variable is 1, the aggregation occurs, 
and on the next MCS the aggregated particles move together as a whole cluster. 

The tick size of the simulation clock at a given MCS, At, is adjusted to the fastest 
particle with Vmax and also scaled by the lattice size and number of dimensions. Since 
the particles start to move faster towards the end of simulation, the clock slows down, 
and it takes more MCS for the latest drying stages. Reasonably, larger lattices also 
require more MCS to finish the whole simulation. The diffusion coefficient, according 
to Codling et al.^^, in BRW model is proportional to sum of individual move prob- 
abilities, 



Px. +Px^ +Py. +Py^ +Pz. +Pz^=l-Pw = 



- (Vx +Vy + Vz) + (Vx + Vy + Vz)Vm, 



(3) 



Therefore, shorter MCS with a larger Vmax have smaller random displacement caused 
by the Brownian motion. The waiting probability, p^, can be alternatively expressed 
as, 



{NVmax - (Vx + Vy + Vz)) {Vmax " 1 ) 



(4) 
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Figure 4 | Analysis of the height profiles of the structures presented in Fig. 3. (al-a5) show the color maps of the particle structure height corresponding 
to particle concentrations, (j), from (al) 0.01 to (a5) 0.05, andpstkh = O-l- (bl-b5) showthe corresponding maps for different sticking p^^cfc = 1- The color 
bar shows the dependence of the height value on the displayed color. Minimal height shown is 0 (empty area), maximum height is 9, meaning that nine 
particle layers are stacked on each other in these cells of the lattice. The plot diagrams (a6) and (b6) show the corresponding ring thickness profiles 
depending on the distance from the domain center. The ring height is averaged over the all cross-sectional directions. His the height of the particle layer, r 
is the non-dimensional distance from the base center, and R is the base radius respectively. 



It shows that for the fastest possible particle with (v^ax, Vmax, Vmax) some move 
definitely happens (p^ = 0); and if all the particles are slow {vmax = l)> there is no 
waiting on the current MCS. In the actual simulations Vmax is almost always higher 
than 1. However, if this condition is not fulfilled, the probabilities are not rescaled by 
Vmax- Then the condition {p^ = 0) is used, and the equations for the slow drift v/vmax 
1 



of the form p = 



are used to calculate pmove > in which V(x,jk,z) is the 



2N - 2NVmax 

projection value on x, y, or z axis respectively. 
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